home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Support / MDPlotting.m < prev    next >
Encoding:
Text File  |  1992-08-18  |  11.7 KB  |  394 lines
  1. (*  :Title:    Supplemental for Plotting Multidimensional Transforms  *)
  2.  
  3. (*  :Authors:    Brian Evans, James McClellan  *)
  4.  
  5. (*
  6.     :Summary:    Plots of m-D magnitude/phase responses
  7.         Generates of m-D pole-zero diagrams and root loci
  8.  *)
  9.  
  10. (*  :Context:    SignalProcessing`Support`TransSupport`  *)
  11.  
  12. (*  :PackageVersion:  2.7    *)
  13.  
  14. (*
  15.     :Copyright:    Copyright 1989-1991 by Brian L. Evans
  16.         Georgia Tech Research Corporation
  17.  
  18.     Permission to use, copy, modify, and distribute this software
  19.     and its documentation for any purpose and without fee is
  20.     hereby granted, provided that the above copyright notice
  21.     appear in all copies and that both that copyright notice and
  22.     this permission notice appear in supporting documentation,
  23.     and that the name of the Georgia Tech Research Corporation,
  24.     Georgia Tech, or Georgia Institute of Technology not be used
  25.     in advertising or publicity pertaining to distribution of the
  26.     software without specific, written prior permission.  Georgia
  27.     Tech makes no representations about the suitability of this
  28.     software for any purpose.  It is provided "as is" without
  29.     express or implied warranty.
  30.  *)
  31.  
  32. (*  :History:    *)
  33.  
  34. (*  :Keywords:    pole-zero diagrams, frequency response  *)
  35.  
  36. (*
  37.     :Source:    {Multidimensional Digital Signal Processing}, 1984,
  38.           by D. Dudgeon and R. Mersereau
  39.  *)
  40.  
  41. (*
  42.     :Warning:    The transform objects LTransData and ZTransData must
  43.         be defined before this file is loaded (see "ROC.m").
  44.  *)
  45.  
  46. (*  :Mathematica Version:  1.2 or 2.0  *)
  47.  
  48. (*  :Limitation:  *)
  49.  
  50. (*
  51.     :Discussion:  This file has two sections:
  52.  
  53.         A.  Magnitude and phase plots (for 2-D signals)
  54.         B.  Pole-zero diagrams and root maps (for 2-D transforms)
  55.  
  56.     These sections display the results of transform objects.
  57.  *)
  58.  
  59. (*  :Functions:      *)
  60.  
  61.  
  62. If [ TrueQ[ $VersionNumber >= 2.0 ],
  63.      Off[ General::spell ];
  64.      Off[ General::spell1 ] ]
  65.  
  66.  
  67. (*  B E G I N     P A C K A G E  *)
  68.  
  69. BeginPackage[ "SignalProcessing`Support`MDPlotting`" ]
  70. EndPackage[]
  71.  
  72. BeginPackage[ "SignalProcessing`Support`TransSupport`",
  73.           "SignalProcessing`Support`ROC`",
  74.           "SignalProcessing`Support`SigProc`",
  75.           "SignalProcessing`Support`SupCode`" ]
  76.  
  77.  
  78. Begin["`Private`"]
  79.  
  80.  
  81. (*  A.  F R E Q U E N C Y     R E S P O N S E  *)
  82.  
  83.  
  84. (*  magPhasePlot  *)
  85.  
  86. fmag[w01_Real, omega1_, w02_Real, omega2_, fresp_, logflag_] :=
  87.     Block [    {value},
  88.         value = Abs[ GetValue[fresp, {omega1, omega2}, {w01, w02}] ];
  89.         If [ logflag, 20 Log[10, value], value ] ]
  90.  
  91. fphase[w01_Real, omega1_, w02_Real, omega2_, fresp_, indegrees_] :=
  92.     Block [    {value},
  93.         value = GetValue[fresp, {omega1, omega2}, {w01, w02}];
  94.         If [ ZeroQ[value], Return[0.] ];
  95.         value = N [ Arg[value] ];
  96.         Chop[ If [ indegrees, value / Degree, value ] ] ]
  97.  
  98.  
  99. (*  Magnitude and phase plots for 2-D signals  *)
  100. magPhasePlot2D[freqresp_, {w1_Symbol, wmin1_, wmax1_},
  101.               {w2_Symbol, wmin2_, wmax2_}, op___] :=
  102.     Block [    {arglist, bogus, degrees, dim1, dim2, fresp, ftemp, logrange,
  103.          magfun, magplot = NullPlot, omega1, omega2, omitplot, options,
  104.          phasefun, phaseplot = NullPlot, plotlabel, plotoptions,
  105.          plotrange, result = 0, rules, w1str, w2str},
  106.  
  107.         (*  Check for invalid arguments  *)
  108.  
  109.         If [ N[wmin1 > wmax1],
  110.              Message[Plot::limits, {w1, wmin1, wmax1}]; Return[Null] ];
  111.         If [ N[wmin2 > wmax2],
  112.              Message[Plot::limits, {w2, wmin2, wmax2}]; Return[Null] ];
  113.  
  114.         (*  Set up for plotting  *)
  115.  
  116.         Off[Power::infy, Infinity::indt];
  117.         w1str = StripPackage[w1];
  118.         w2str = StripPackage[w2];
  119.         options = ToList[op] ~Join~ Options[MagPhasePlot];
  120.         plotoptions = RemoveOptions[options, badoptions[MagPhasePlot]];
  121.  
  122.         (*  Adjust frequency response so that plot shows periodicity *)
  123.  
  124.         fresp = freqresp /. PeriodicRules;
  125.  
  126.                 (*  Extract delta functions from the freq. resp.  *)
  127.  
  128.         deltafuns = GetDeltaFunctions[funct, {w1, w2}];
  129.         If [ ! SameQ[deltafuns, 0],
  130.              Message[ MagPhasePlot::nodeltas ];
  131.              fresp = fresp //.
  132.                   { h_[ c_. Delta[a_. w1 + b_.] ] :> 0,
  133.                 h_[ c_. Delta[a_. w1 + b_.] + x_ ] :> h[x],
  134.                 Delta[a_. w1 + b_.] :> 0,
  135.                 h_[ c_. Delta[a_. w2 + b_.] ] :> 0,
  136.                 h_[ c_. Delta[a_. w2 + b_.] + x_ ] :> h[x],
  137.                 Delta[a_. w2 + b_.] :> 0 } ];
  138.  
  139.                 (*  Convert expression to plottable/functional form    *)
  140.         (*  and replace constants like Pi with numerical value *)
  141.  
  142.         fresp = N [ TheFunction[fresp] ];
  143.  
  144.         (*  Find valueless symbols other than frequency variables *)
  145.                 (*  The first parameter in Summation is a local symbol   *)
  146.  
  147.         ftemp = fresp /. 
  148.                           ( Summation[n_, l_, u_, inc_][expr_] :>
  149.                 bogus[l, u, inc, GetVariables[expr /. n -> l]] );
  150.         varlist = Select[GetVariables[ftemp],
  151.                  ((! SameQ[#1, w1]) && (! SameQ[#1, w2]))&];
  152.         If [ ! EmptyQ[varlist],
  153.              Message[ MagPhasePlot::unresolved, varlist];
  154.              fresp = fresp /. Map[(#1 -> 1)&, varlist] ];
  155.  
  156.         (*  Substitute dummy variables  *)
  157.  
  158.         fresp = fresp /. { w1 -> omega1, w2 -> omega2 };
  159.  
  160.         (*  Magnitude plotting  *)
  161.  
  162.         omitplot = SameQ[Replace[MagRangeScale, options], Null];
  163.         logrange = SameQ[Replace[MagRangeScale, options], Log];
  164.         plotlabel = If [ logrange,
  165.                  "Magnitude Response (dB)",
  166.                  "Magnitude Response" ];
  167.         plotrange = Replace[PlotRange, options];
  168.  
  169.         If [ ! omitplot,
  170.              arglist = { fmag[w1, omega1, w2, omega2, fresp, logrange],
  171.                  { w1, wmin1, wmax1 },
  172.                  { w2, wmin2, wmax2 } } ~Join~
  173.                    plotoptions ~Join~
  174.                    { AxesLabel -> { w1str, w2str, " " },
  175.                  DisplayFunction -> disp,
  176.                  PlotLabel -> plotlabel,
  177.                  PlotRange -> plotrange };
  178.              magplot = Apply[ Plot3D, arglist ] ];
  179.  
  180.         (*  Phase plotting  *)
  181.  
  182.         omitplot = SameQ[Replace[PhaseRangeScale, options], Null];
  183.         degrees = SameQ[Replace[PhaseRangeScale, options], Degree];
  184.         plotlabel = If [ degrees,
  185.                  "Phase Response (degrees)",
  186.                  "Phase Response (radians)" ];
  187.  
  188.         If [ ! omitplot,
  189.              arglist = { fphase[w1, omega1, w2, omega2, fresp, degrees],
  190.                  { w1, wmin1, wmax1 },
  191.                  { w2, wmin2, wmax2 } } ~Join~
  192.                    plotoptions ~Join~
  193.                    { DisplayFunction -> disp,
  194.                  PlotLabel -> plotlabel,
  195.                  PlotRange -> All,
  196.                  AxesLabel -> { w1str, w2str, " " } };
  197.              phaseplot = Apply[ Plot3D, arglist ] ];
  198.  
  199.         (*  Clean up and return the plots as graphics objects  *)
  200.  
  201.         On[Power::infy, Infinity::indt];
  202.         { magplot, phaseplot } ]
  203.  
  204.  
  205. (*  B.  P O L E - Z E R O     P L O T T I N G  *)
  206.  
  207. (*      1.  D r i v e r s  f o r  b o t h  z -  a n d  s - d o m a i n s  *)
  208.  
  209. (*  Two-dimensional pole-zero plotting driver                 *)
  210. (*  --  separable transforms generate two separable pole-zero plots  *) 
  211. (*  --  symmetric transfer function f only requires one plot         *)
  212. (*  --  otherwise, project z1 onto z2 and z2 onto z1             *)
  213. poleZeroPlot2D[f_, {z1_Symbol, z2_Symbol},
  214.            {rm1_, rm2_}, {rp1_, rp2_}, rest__ ] :=
  215.     Block [    {seplist, z1expr, z2expr},
  216.         If [ SeparableQ[f, {z1, z2}],
  217.              Block [ {poles, seplist, z1expr, z2expr},
  218.                  seplist = Separate[f, {z1, z2}];
  219.                  z1expr = seplist[[1]];
  220.                  z2expr = seplist[[2]];
  221.                  poles = PoleZeroPlot[z1expr, z1, rm1, rp1, rest];
  222.                  poles = poles ~Join~
  223.                      PoleZeroPlot[z2expr, z2, rm2, rp2, rest];
  224.                  poles ],
  225.              If [ ! SameQ[f, f /. { z1 -> z2, z2 -> z1 }],
  226.                   poleZeroRootMap[f, z2, z1, rm1, rp1, rest] ];
  227.              poleZeroRootMap[f, z1, z2, rm2, rp2, rest] ] ]
  228.  
  229. PoleZeroRootLocus[ ztrans_, z_Symbol, {w_Symbol, wmin_, wmax_, wstep_},
  230.            ucROCplot_, options_ ] :=
  231.     Block [    {denom, numer, polelist, rmapops, rootmap,
  232.          showops, xlabel, ylabel, zerolist, zstr},
  233.  
  234.         numer = Numerator[ztrans];
  235.         denom = Denominator[ztrans];
  236.  
  237.         (*  Tell user what we're doing if Dialogue is enabled  *)
  238.  
  239.         If [ InformUserQ[options],
  240.              Print[ " " ];
  241.              Print[ "Numerator polynomial in ", z, ": ",
  242.                 numer /. w -> "w" ];
  243.              Print[ "Denominator polynomial in ", z, ": ",
  244.                 denom /. w -> "w" ] ];
  245.  
  246.         (*  Root map labels and options  *)
  247.  
  248.         zstr = StripPackage[z];
  249.         xlabel = "Re " ~StringJoin~ zstr;
  250.         ylabel = "Im " ~StringJoin~ zstr;
  251.         showops = Append[ RemoveOptions[options, {Dialogue}],
  252.                   AxesLabel -> { xlabel, ylabel } ];
  253.         rmapops = Prepend[ showops, DisplayFunction :> Identity ];
  254.  
  255.         (*  Zero root map  *)
  256.  
  257.         If [ ! FreeQ[numer, z],
  258.              rootmap = RootLocus[ numer, z, {w, wmin, wmax, wstep},
  259.                            rmapops ][[2]];
  260.              Show[ ucROCplot, rootmap,
  261.                Append[showops, PlotLabel -> "Zero Root Map"] ] ];
  262.  
  263.         (*  Pole root map  *)
  264.  
  265.         If [ ! FreeQ[denom, z],
  266.              rootmap = RootLocus[ 1/denom, z, {w, wmin, wmax, wstep},
  267.                       rmapops ][[2]];
  268.              Show[ ucROCplot, rootmap,
  269.                Append[showops, PlotLabel -> "Pole Root Map"] ] ];
  270.  
  271.         GetRootList[denom, z] ]
  272.  
  273.  
  274. (*      2.  P l o t t i n g   r o u t i n e s   f o r   z - d o m a i n  *)
  275.  
  276. (*  Root map two-dimensional pole-zero plotting    *)
  277. (*  variable zproj is projected onto z-domain    *)
  278. poleZeroRootMap[f_, zproj_Symbol, z_Symbol, rm_, rp_, True, op___] :=
  279.     Block [ {graphfuns, options, plotstyle,
  280.          rminus, rplus, ucROCplot, w, ztrans},
  281.  
  282.         options = pzOptions[op] ~Join~
  283.               { AspectRatio -> 1, PlotRange -> {{-2, 2}, {-2, 2}} };
  284.  
  285.         (*  Make sure that the transform is over common denominator  *)
  286.  
  287.         ztrans = Together[f];
  288.         If [ ! MyRationalPolyQ[ztrans, z],
  289.              Message[ PoleZeroPlot::notrational ];
  290.              Message[ PoleZeroPlot::noplot ];
  291.              Return[{}] ];
  292.  
  293.         (*  We will plot a unit circle and R- and possibly R+  *)
  294.  
  295.         rminus = N[rm /. zproj -> w];
  296.         rplus = N[rp /. zproj -> w];
  297.  
  298.         plotstyle = { Thickness[0.01] };
  299.         graphfuns = { { Cos[w], Sin[w] } };
  300.  
  301.         If [ ! ZeroQ[rminus],
  302.              PrependTo[ graphfuns, { rminus Cos[w], rminus Sin[w] } ];
  303.              PrependTo[ plotstyle, Dashing[{0.05,0.05}] ] ];
  304.  
  305.         If [ ! InfinityQ[rplus],
  306.              PrependTo[ graphfuns, { rplus Cos[w], rplus Sin[w] } ];
  307.              PrependTo[ plotstyle, Dashing[{0.05,0.05}] ] ];
  308.  
  309.         ucROCplot = Apply[ ParametricPlot,
  310.                    { graphfuns,
  311.                      {w, 0, 2 Pi},
  312.                      AspectRatio -> 1,
  313.                      DisplayFunction :> Identity,
  314.                      PlotStyle -> plotstyle } ];
  315.  
  316.         (*  Project the transfer function and draw the root locus  *)
  317.  
  318.         ztrans = ztrans /. zproj -> Exp[I w];
  319.  
  320.         PoleZeroRootLocus[ ztrans, z, {w, 0, 2 Pi, 2 Pi / 20},
  321.                    ucROCplot, options ] ]
  322.  
  323.  
  324. (*      3.  P l o t t i n g   R o u t i n e s   f o r   s - d o m a i n  *)
  325.  
  326. (*  Root map two-dimensional pole-zero plotting    *)
  327. (*  variable sproj is projected onto s-domain    *)
  328. poleZeroRootMap[f_, sproj_Symbol, s_Symbol, rm_, rp_, False, op___] :=
  329.     Block [ {graphfuns, ltrans, options, plotstyle, rminus,
  330.          rplus, ucROCplot, ucROCplot1, ucROCplot2, w},
  331.  
  332.         options = pzOptions[op];
  333.  
  334.         (*  Make sure the transfer function is a ratio of two polys  *)
  335.  
  336.         ltrans = Together[f];
  337.         If [ ! MyRationalPolyQ[ltrans, s],
  338.              Message[ PoleZeroPlot::notrational ];
  339.              Message[ PoleZeroPlot::noplot ];
  340.              Return[{}] ];
  341.  
  342.         (*  We will possibly plot R- and R+  *)
  343.  
  344.         rminusline = NullPlot;
  345.         rplusline = NullPlot;
  346.  
  347.         rminus = N[rm /. sproj -> (I w)];
  348.         rplus = N[rp /. sproj -> (I w)];
  349.  
  350.         plotstyle = { };
  351.         graphfuns = { };
  352.  
  353.         If [ ! InfinityQ[rminus],
  354.              PrependTo[ graphfuns, {rminus, w} ];
  355.              PrependTo[ plotstyle, Dashing[{0.05,0.05}] ] ];
  356.  
  357.         If [ ! InfinityQ[rplus],
  358.              PrependTo[ graphfuns, {rplus, w} ];
  359.              PrependTo[ plotstyle, Dashing[{0.05,0.05}] ] ];
  360.  
  361.         ucROCplot1 = Apply[ ParametricPlot,
  362.                     { graphfuns,
  363.                       {w, -10, 10},
  364.                       DisplayFunction :> Identity,
  365.                       PlotStyle -> plotstyle } ];
  366.  
  367.         ucROCplot2 = Apply[ ParametricPlot,
  368.                     { graphfuns,
  369.                       {w, -10000, 10000},
  370.                       DisplayFunction :> Identity,
  371.                       PlotStyle -> plotstyle } ];
  372.  
  373.         ucROCplot = Show[ {ucROCplot1, ucROCplot2},
  374.                   DisplayFunction :> Identity ];
  375.  
  376.         (*  Project the transfer function and draw the root locus  *)
  377.  
  378.         ltrans = ltrans /. sproj -> (I w);
  379.  
  380.         PoleZeroRootLocus[ ltrans, s, {w, -10000, 10000, 1000},
  381.                    ucROCplot, options ] ]
  382.  
  383.  
  384. (*  E N D     P A C K A G E  *)
  385.  
  386. End[]
  387. EndPackage[]
  388.  
  389.  
  390. (*  E N D I N G     M E S S A G E  *)
  391.  
  392. Print["Multidimensional plotting routines for transforms are loaded."]
  393. Null
  394.